Results for Hitchcock & Frank 2023
Peter Frank Hitchcock - LNCC Lab, Brown University
2023-06-23
- Load data from studies 1 and 2
- Behavioral Results for Learning and Test phase
- Modeling Results (Note: model numbers match those from
s.Rbut not those used in paper…)- Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter
- Model 3 — Same as m1 but not letting decay vary
- Model 4 — Same as m3 but to 0 inits
- Model 5 — Same as m1 but lets beta vary instead of decay
- Model 6 — Same as m1 but lets learning rate vary instead of decay
- Model 11 — Simple q-learner model
- Model 12 — Simple Bayes model
- Model 13 — M12 but allowing beta to vary
- Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond
- Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond
- Model 27 — Fixing both ES and epsilon at the group level
- Model 29 — Model 1 but with no choice kernel
- Model Comparison
sapply(c(
"rjson",
"data.table",
"dplyr",
"ggplot2",
"stringr",
"purrr",
"foreach",
"doParallel",
"patchwork",
"lme4",
"lmerTest",
"testit",
"ggpubr",
"latex2exp"
),
require, character=TRUE)## Loading required package: rjson## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: ggplot2## Loading required package: stringr## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: foreach##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, when## Loading required package: doParallel## Loading required package: iterators## Loading required package: parallel## Loading required package: patchwork## Loading required package: lme4## Loading required package: Matrix## Loading required package: lmerTest##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step## Loading required package: testit## Loading required package: ggpubr## Loading required package: latex2exp## rjson data.table dplyr ggplot2 stringr purrr foreach
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## doParallel patchwork lme4 lmerTest testit ggpubr latex2exp
## TRUE TRUE TRUE TRUE TRUE TRUE TRUEsf <- function() sapply(paste0("./Functions/", list.files("./Functions/", recursive=TRUE)), source) # Source all fxs
sf()## ./Functions/General/FindDelay.R ./Functions/General/FindDelayAlt.R
## value ? ?
## visible FALSE FALSE
## ./Functions/General/Utils.R ./Functions/Model/Models/RunMBayesLearner.R
## value ? ?
## visible FALSE FALSE
## ./Functions/Model/Models/RunMBayesLearnerDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDecayingQLearnerDiffBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayes.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMHybridDiffDecayQLearnerBayesAlt.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearner.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayTo0Inits.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffBeta.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDecayToPessPriorDiffLR.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPrior.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorESAndEpsFixed.R
## value ?
## visible FALSE
## ./Functions/Model/Models/RunMQLearnerDiffDecayToPessPriorNoCK.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/AltPlotSimTrainingCurves.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/aModelUtils.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotAllWorstTestOptionsSimVsEmp.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTest.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PlotSimEmpTestNoSim.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModel.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/PrepDfForModelPROpt.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/RecodeDfs.R
## value ?
## visible FALSE
## ./Functions/Model/ModelUtils/SimplePlotSimTrainingDelay.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/OptFxs/RunOptTrainSITPR.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSIT.R
## value ?
## visible FALSE
## ./Functions/Model/SimFxs/RunSimTrainSITForPR.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/aModelFunctions.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesMean.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesStd.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/CalcBayesVar.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayChoiceKernel.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/DecayQVals.R
## value ?
## visible FALSE
## ./Functions/Model/TrialWiseComps/UpdateABMatrix.R
## value ?
## visible FALSEDefPlotPars()
registerDoParallel(cores=round(detectCores()*2/3))Load data from studies 1 and 2
s1_train <- read.csv("../data/cleaned_files/s1_train_with_delay_deident.csv")
#read.csv("../data/cleaned_files/s1_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s1_sit <- read.csv("../data/cleaned_files/s1_SIT_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_train <- read.csv("../data/cleaned_files/s2_train_with_delay_deident.csv") #read.csv("../data/cleaned_files/s2_train_deident.csv") %>% rename(ID=deident_ID) %>% relocate(ID)
s2_sit <- read.csv("../data/cleaned_files/s2_sit_deident_corrected_names.csv") Harmonize variables and create some separate vars
# Study 2 harmonize
s2_sit[s2_sit$condition=="cogn", "condition"] <- "cognitive"
s2_train[s2_train$trial_within_condition <= 20, "block"] <- 1
s2_train[s2_train$trial_within_condition > 20, "block"] <- 2
s1_sit$probability <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 1)))
s1_sit$valence <- factor(unlist(map(strsplit(as.character(s1_sit$valence_and_probability), "_"), 2)))
s2_sit$probability <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 1)))
s2_sit$valence <- factor(unlist(map(strsplit(as.character(s2_sit$valence_and_probability), "_"), 2)))Define paths and functions
# MLE results
allp <- "../model_res/opts_mle_paper_final/all/"
# Empirical Bayes results
bp <- "../model_res/opts_mle_paper_final/best/"
# Simulations
sp <- "../model_res/sims_clean/sims_from_mle/"
# Read model function
rm <- function(path, model_str) read.csv(paste0(path, model_str))Behavioral Results for Learning and Test phase
Visualization
Summarize training mean and within-subject SEM
s1_train_summs <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_tr_summs <- Rmisc::summarySEwithin(s1_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_train_summs <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_tr_summs <- Rmisc::summarySEwithin(s2_train_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions1_test_summs <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_te_summs <- Rmisc::summarySEwithin(s1_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_test_summs <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions2_te_summs <- Rmisc::summarySEwithin(s2_test_summs,
measurevar = "m",
withinvars = "condition",
idvar = "ID")## Automatically converting the following non-factors to factors: conditions1_tr_summs$experiment <- 1
s1_tr_summs$phase <- "Learning"
s2_tr_summs$experiment <- 2
s2_tr_summs$phase <- "Learning"
s1_te_summs$experiment <- 1
s1_te_summs$phase <- "Test"
s2_te_summs$experiment <- 2
s2_te_summs$phase <- "Test"
comb_summs <- rbind(s1_tr_summs, s2_tr_summs)
comb_summs$experiment <- factor(comb_summs$experiment)
comb_te_summs <- rbind(s1_te_summs, s2_te_summs)
comb_te_summs$experiment <- factor(comb_te_summs$experiment)a <- ggplot(comb_summs,
aes(x=experiment, y=m, group=condition, fill=condition)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) +
geom_hline(yintercept=.5, size=2) +
scale_fill_manual(values=c("purple", "orange")) +
geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
group=condition),
inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
geom_point(size=4, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol +
xlab("") + ylab("Proportion correct") + ggtitle("Learning") + tp +
ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.a
b <- ggplot(comb_te_summs,
aes(x=experiment, y=m, group=condition, fill=condition)) +
geom_hline(yintercept=c(seq(.6, .8, .1)), size=.3) +
geom_hline(yintercept=.5, size=2) +
scale_fill_manual(values=c("purple", "orange")) +
geom_errorbar(aes(x=experiment, ymin=m-se, ymax=m+se,
group=condition),
inherit.aes=FALSE, size=1.5, width=.25, position = position_dodge(width = .3)) +
geom_point(size=3, color="black", pch=21, position = position_dodge(width = .3)) + ga + ap + tol + tp +
#lp + Used to cut legend and then turned off so size is comparable
xlab("") + ylab("") + ggtitle("Test") +
ylim(c(.45, .82)) + theme(axis.text.x = element_text(size=30))
b
perf_comb <- a + bperf_comb
#ggsave("../paper/figs/fig_parts/perf_plot.png", perf_comb, height=6, width=10, dpi=700)Statistical analyses
Create sum contrast codes
# Specify overt as baseline factor so that negative effect corresponds to worst performance
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
#head(s1_train$cond_cc)
s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)
s1_sit$cond_cc <- factor(s1_sit$condition, levels=c("overt", "cognitive"))
contrasts(s1_sit$cond_cc) <- c(-.5, .5)
s2_sit$cond_cc <- factor(s2_sit$condition, levels=c("overt", "cognitive"))
contrasts(s2_sit$cond_cc) <- c(-.5, .5)Regressions of condition differences in performance
Experiment 1 — Learning and test models
summary(m1_train_full_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +
## scale(trial_within_condition) | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45054.2 45131.6 -22518.1 45036.2 39991
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.9157 -1.0686 0.4573 0.7047 1.2827
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4466 0.6683
## cond_cc1 0.4926 0.7018 -0.19
## scale(trial_within_condition) 0.1156 0.3399 0.82 0.01
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.04191 0.06122 17.019 < 2e-16 ***
## cond_cc1 -0.18054 0.06745 -2.677 0.00743 **
## scale(trial_within_condition) 0.36991 0.03300 11.210 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.174
## scl(trl_w_) 0.758 0.009summary(m1_test_full_regr <- glmer(correct ~ cond_cc + (cond_cc|ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7854.4 7889.3 -3922.2 7844.4 7995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1729 0.1319 0.2789 0.5791 1.2316
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.558 1.248
## cond_cc1 1.972 1.404 -0.17
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.5690 0.1193 13.149 < 2e-16 ***
## cond_cc1 -0.4569 0.1510 -3.027 0.00247 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.173Experiment 2 — Learning and test models
summary(m2_train_full_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(cond_cc + scale(trial_within_condition) | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (cond_cc +
## scale(trial_within_condition) | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 48366.7 48445.0 -24174.4 48348.7 44151
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.0445 -1.0336 0.4405 0.6872 1.4792
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6677 0.8171
## cond_cc1 0.5648 0.7515 -0.07
## scale(trial_within_condition) 0.1676 0.4094 0.87 -0.04
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.14913 0.07101 16.184 <2e-16 ***
## cond_cc1 -0.17464 0.06852 -2.549 0.0108 *
## scale(trial_within_condition) 0.39251 0.03733 10.515 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.071
## scl(trl_w_) 0.823 -0.034summary(m2_test_full_regr <- glmer(correct ~ cond_cc + (cond_cc|ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8426.3 8461.7 -4208.1 8416.3 8827
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2283 0.1250 0.2301 0.5751 1.5541
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.850 1.360
## cond_cc1 2.467 1.571 -0.15
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.6720 0.1235 13.542 <2e-16 ***
## cond_cc1 -0.3543 0.1586 -2.234 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.156car::vif(m1_train_full_regr)## cond_cc scale(trial_within_condition)
## 1.000075 1.000075car::vif(m2_train_full_regr)## cond_cc scale(trial_within_condition)
## 1.001187 1.001187Robustness check 1: Same models with just random intercepts instead of both random intercepts and slopes
Experiment 1
summary(m1_train_ri_only_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(1 |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46071.1 46105.5 -23031.6 46063.1 39996
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6755 -1.1011 0.5027 0.6818 1.2499
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.322 0.5675
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95856 0.05214 18.386 < 2e-16 ***
## cond_cc1 -0.14312 0.02271 -6.302 2.93e-10 ***
## scale(trial_within_condition) 0.28711 0.01147 25.021 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.006
## scl(trl_w_) 0.026 -0.008summary(m1_test_ri_only_regr <- glmer(correct ~ cond_cc + (1 |ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8080.6 8101.5 -4037.3 8074.6 7997
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4964 0.1516 0.3789 0.6174 1.2182
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.299 1.14
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.40802 0.10812 13.023 < 2e-16 ***
## cond_cc1 -0.31767 0.05567 -5.707 1.15e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.022Experiment 2
summary(m2_train_ri_only_regr <- glmer(correct ~
cond_cc + scale(trial_within_condition) +
(1 | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + scale(trial_within_condition) + (1 | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 49635.6 49670.3 -24813.8 49627.6 44156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4739 -1.0679 0.4853 0.6668 1.2916
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.4132 0.6428
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.03285 0.05596 18.458 < 2e-16 ***
## cond_cc1 -0.15224 0.02195 -6.937 4.01e-12 ***
## scale(trial_within_condition) 0.27618 0.01108 24.922 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1
## cond_cc1 -0.007
## scl(trl_w_) 0.024 -0.008summary(m2_test_ri_only_regr <- glmer(correct ~ cond_cc + (1 |ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (1 | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8769.8 8791.0 -4381.9 8763.8 8829
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6214 0.1441 0.3650 0.6075 1.2670
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 1.577 1.256
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.48452 0.11314 13.121 < 2e-16 ***
## cond_cc1 -0.20840 0.05335 -3.906 9.38e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.013Robustness check 2: Learning models with the time covariate removed
Experiments 1 and 2
summary(m1_train_no_cov_regr <- glmer(correct ~
cond_cc +
(cond_cc |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46182.9 46225.9 -23086.5 46172.9 39995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3753 -1.1134 0.5036 0.6746 1.0940
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.3342 0.5781
## cond_cc1 0.4616 0.6794 -0.20
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.97023 0.05303 18.296 < 2e-16 ***
## cond_cc1 -0.17611 0.06533 -2.696 0.00702 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.189summary(m2_train_no_cov_regr <- glmer(correct ~
cond_cc +
(cond_cc | ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ cond_cc + (cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 49610.0 49653.5 -24800.0 49600.0 44155
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9719 -1.0789 0.4729 0.6802 1.0720
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.4287 0.6548
## cond_cc1 0.5173 0.7192 -0.09
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.04886 0.05698 18.407 < 2e-16 ***
## cond_cc1 -0.16968 0.06571 -2.582 0.00981 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cond_cc1 -0.088sjPlot::tab_model(m1_train_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.83 | 2.51 – 3.20 | <0.001 |
| cond cc1 | 0.83 | 0.73 – 0.95 | 0.007 |
| trial within condition | 1.45 | 1.36 – 1.54 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.45 | ||
| τ11 ID.cond_cc1 | 0.49 | ||
| τ11 ID.scale(trial_within_condition) | 0.12 | ||
| ρ01 | -0.19 | ||
| 0.82 | |||
| ICC | 0.17 | ||
| N ID | 125 | ||
| Observations | 40000 | ||
| Marginal R2 / Conditional R2 | 0.035 / 0.202 | ||
sjPlot::tab_model(m2_train_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 3.16 | 2.75 – 3.63 | <0.001 |
| cond cc1 | 0.84 | 0.73 – 0.96 | 0.011 |
| trial within condition | 1.48 | 1.38 – 1.59 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 0.67 | ||
| τ11 ID.cond_cc1 | 0.56 | ||
| τ11 ID.scale(trial_within_condition) | 0.17 | ||
| ρ01 | -0.07 | ||
| 0.87 | |||
| ICC | 0.23 | ||
| N ID | 138 | ||
| Observations | 44160 | ||
| Marginal R2 / Conditional R2 | 0.037 / 0.257 | ||
sjPlot::tab_model(m1_test_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 4.80 | 3.80 – 6.07 | <0.001 |
| cond cc1 | 0.63 | 0.47 – 0.85 | 0.002 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.56 | ||
| τ11 ID.cond_cc1 | 1.97 | ||
| ρ01 ID | -0.17 | ||
| ICC | 0.38 | ||
| N ID | 125 | ||
| Observations | 8000 | ||
| Marginal R2 / Conditional R2 | 0.010 / 0.390 | ||
sjPlot::tab_model(m2_test_full_regr)| correct | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 5.32 | 4.18 – 6.78 | <0.001 |
| cond cc1 | 0.70 | 0.51 – 0.96 | 0.025 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 1.85 | ||
| τ11 ID.cond_cc1 | 2.47 | ||
| ρ01 ID | -0.15 | ||
| ICC | 0.43 | ||
| N ID | 138 | ||
| Observations | 8832 | ||
| Marginal R2 / Conditional R2 | 0.005 / 0.432 | ||
Robustness check 3: Covariate for answer type
s1_train[s1_train$response=="sum", "resp_type_cc_cog"] <- .5
s1_train[s1_train$response=="difference", "resp_type_cc_cog"] <- -.5
s1_train[s1_train$response=="top", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="bottom", "resp_type_cc_cog"] <- 0
s1_train[s1_train$response=="top", "resp_type_cc_overt"] <- .5
s1_train[s1_train$response=="bottom", "resp_type_cc_overt"] <- -.5
s1_train[s1_train$response=="sum", "resp_type_cc_overt"] <- 0
s1_train[s1_train$response=="difference", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_cog"] <- .5
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_cog"] <- -.5
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_cog"] <- 0
s1_sit[s1_sit$resp_as_category=="top", "resp_type_cc_overt"] <- .5
s1_sit[s1_sit$resp_as_category=="bottom", "resp_type_cc_overt"] <- -.5
s1_sit[s1_sit$resp_as_category=="sum", "resp_type_cc_overt"] <- 0
s1_sit[s1_sit$resp_as_category=="difference", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="alphabetize", "resp_type_cc_cog"] <- .5
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_train[s2_train$response=="slash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="backslash", "resp_type_cc_cog"] <- 0
s2_train[s2_train$response=="slash", "resp_type_cc_overt"] <- .5
s2_train[s2_train$response=="backslash", "resp_type_cc_overt"] <- -.5
s2_train[s2_train$response=="alphabetize", "resp_type_cc_overt"] <- 0
s2_train[s2_train$response=="rev_alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_cog"] <- .5
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_cog"] <- -.5
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_cog"] <- 0
s2_sit[s2_sit$resp_as_category=="slash", "resp_type_cc_overt"] <- .5
s2_sit[s2_sit$resp_as_category=="backslash", "resp_type_cc_overt"] <- -.5
s2_sit[s2_sit$resp_as_category=="alphabetize", "resp_type_cc_overt"] <- 0
s2_sit[s2_sit$resp_as_category=="rev_alphabetize", "resp_type_cc_overt"] <- 0
table(s1_train$resp_type_cc_cog)##
## -0.5 0 0.5
## 9020 20000 10980table(s1_train$resp_type_cc_overt)##
## -0.5 0 0.5
## 9451 20000 10549table(s1_sit$resp_type_cc_cog)##
## -0.5 0 0.5
## 1851 4000 2149table(s1_sit$resp_type_cc_overt)##
## -0.5 0 0.5
## 1975 4000 2025table(s2_train$resp_type_cc_cog)##
## -0.5 0 0.5
## 10182 22080 11898table(s2_train$resp_type_cc_overt)##
## -0.5 0 0.5
## 11186 22080 10894table(s2_sit$resp_type_cc_cog)##
## -0.5 0 0.5
## 2024 4416 2392table(s2_sit$resp_type_cc_overt)##
## -0.5 0 0.5
## 2239 4416 2177Separately coded contrast variables for cog and overt. Random slopes of these removed for test bc singular (which makes sense bc of few observations/subj)
summary(m1_train_full_regr_at <- glmer(correct ~
cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +
## resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +
## resp_type_cc_cog + resp_type_cc_overt | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45040 45212 -22500 45000 39980
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.2191 -1.0643 0.4553 0.7031 1.4380
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.44641 0.6681
## cond_cc1 0.49626 0.7045 -0.19
## scale(trial_within_condition) 0.11598 0.3406 0.81 0.01
## resp_type_cc_cog 0.05454 0.2335 0.19 -0.31 0.33
## resp_type_cc_overt 0.07933 0.2817 -0.04 0.23 0.16 0.70
## Number of obs: 40000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.06182 0.06141 17.290 < 2e-16 ***
## cond_cc1 -0.18291 0.06839 -2.674 0.00749 **
## scale(trial_within_condition) 0.37205 0.03308 11.248 < 2e-16 ***
## resp_type_cc_cog -0.07227 0.04151 -1.741 0.08166 .
## resp_type_cc_overt -0.08409 0.04532 -1.855 0.06353 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1 -0.172
## scl(trl_w_) 0.753 0.004
## rsp_typ_cc_c 0.077 -0.175 0.164
## rsp_typ_cc_v -0.034 0.137 0.088 0.197# Singular
# summary(m1_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
# (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID),
# data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m1_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc |ID),
data=s1_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc | ID)
## Data: s1_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 7857.3 7906.2 -3921.6 7843.3 7993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0858 0.1296 0.2776 0.5823 1.2399
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.555 1.247
## cond_cc1 1.979 1.407 -0.17
## Number of obs: 8000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.57061 0.11923 13.173 < 2e-16 ***
## cond_cc1 -0.45590 0.15120 -3.015 0.00257 **
## resp_type_cc_cog -0.05969 0.08348 -0.715 0.47460
## resp_type_cc_overt -0.07056 0.08936 -0.790 0.42974
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1 -0.175
## rsp_typ_cc_c -0.015 -0.021
## rsp_typ_cc_v -0.007 0.011 0.000summary(m2_train_full_regr_at <- glmer(correct ~
cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc + scale(trial_within_condition) + resp_type_cc_cog + resp_type_cc_overt |ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + scale(trial_within_condition) + resp_type_cc_cog +
## resp_type_cc_overt + (cond_cc + scale(trial_within_condition) +
## resp_type_cc_cog + resp_type_cc_overt | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 48297.6 48471.5 -24128.8 48257.6 44140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.4618 -1.0313 0.4323 0.6859 1.4685
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6757 0.8220
## cond_cc1 0.5985 0.7736 -0.07
## scale(trial_within_condition) 0.1666 0.4081 0.87 -0.04
## resp_type_cc_cog 0.2001 0.4474 0.01 -0.26 0.01
## resp_type_cc_overt 0.1730 0.4160 -0.27 0.00 -0.27 0.18
## Number of obs: 44160, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.19401 0.07162 16.671 <2e-16 ***
## cond_cc1 -0.16320 0.07110 -2.296 0.0217 *
## scale(trial_within_condition) 0.39842 0.03726 10.694 <2e-16 ***
## resp_type_cc_cog -0.12226 0.05288 -2.312 0.0208 *
## resp_type_cc_overt -0.00198 0.05164 -0.038 0.9694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 sc(__) rsp_typ_cc_c
## cond_cc1 -0.064
## scl(trl_w_) 0.821 -0.037
## rsp_typ_cc_c -0.001 -0.199 0.016
## rsp_typ_cc_v -0.183 0.007 -0.170 0.087# Singular
# summary(m2_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
# (cond_cc + resp_type_cc_cog + resp_type_cc_overt |ID),
# data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m2_test_full_regr_at <- glmer(correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt +
(cond_cc |ID),
data=s2_sit, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## correct ~ cond_cc + resp_type_cc_cog + resp_type_cc_overt + (cond_cc | ID)
## Data: s2_sit
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 8427.9 8477.5 -4207.0 8413.9 8825
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3845 0.1211 0.2303 0.5752 1.5798
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 1.848 1.359
## cond_cc1 2.469 1.571 -0.16
## Number of obs: 8832, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.67413 0.12343 13.563 <2e-16 ***
## cond_cc1 -0.35472 0.15864 -2.236 0.0254 *
## resp_type_cc_cog -0.05734 0.08376 -0.685 0.4936
## resp_type_cc_overt 0.11865 0.08580 1.383 0.1667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cnd_c1 rsp_typ_cc_c
## cond_cc1 -0.159
## rsp_typ_cc_c -0.014 -0.022
## rsp_typ_cc_v 0.010 -0.015 0.000Robustness check 4: paired t-test on arcsine-transformed proportion data instead of mixed-effects model
s1_summs_id <- s1_train %>% group_by(condition, ID) %>%
summarize(m=mean(correct)) ## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_summs_id <- s2_train %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s1_summs_test_id <- s1_sit %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.s2_summs_test_id <- s2_sit %>% group_by(condition, ID) %>%
summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.t.test(asin(sqrt(m)) ~ condition, data=s1_summs_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.5152, df = 124, p-value = 0.01318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.061942537 -0.007385729
## sample estimates:
## mean of the differences
## -0.03466413t.test(asin(sqrt(m)) ~ condition, data=s2_summs_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.5139, df = 137, p-value = 0.0131
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.061281292 -0.007319414
## sample estimates:
## mean of the differences
## -0.03430035t.test(asin(sqrt(m)) ~ condition, data=s1_summs_test_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -2.6957, df = 124, p-value = 0.008001
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.12657517 -0.01939599
## sample estimates:
## mean of the differences
## -0.07298558t.test(asin(sqrt(m)) ~ condition, data=s2_summs_test_id, paired=TRUE)##
## Paired t-test
##
## data: asin(sqrt(m)) by condition
## t = -1.8641, df = 137, p-value = 0.06445
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.105733994 0.003119297
## sample estimates:
## mean of the differences
## -0.05130735Delay effects
Create sum contrast codes for delay
# Specify overt as baseline factor so that negative effect corresponds to worst performance
s1_train$cond_cc <- factor(s1_train$condition, levels=c("overt", "cognitive"))
contrasts(s1_train$cond_cc) <- c(-.5, .5)
s2_train$cond_cc <- factor(s2_train$condition, levels=c("overt", "cognitive"))
contrasts(s2_train$cond_cc) <- c(-.5, .5)#head(factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay")))
s1_train$delay_cc <- factor(if_else(s1_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s1_train$delay_cc) <- c(-.5, .5)
s2_train$delay_cc <- factor(if_else(s2_train$delay==0, "no_delay", "delay"), levels=c("no_delay", "delay"))
contrasts(s2_train$delay_cc) <- c(-.5, .5)Building up in complexity…
Effect of delay with no moderation
summary(m1_delay <- glmer(correct ~ delay_cc + (delay_cc|ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 43647.9 43690.6 -21818.9 43637.9 37995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2305 -1.1242 0.5124 0.6818 1.0042
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.32784 0.5726
## delay_cc1 0.02447 0.1564 0.43
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.07222 0.05325 20.14 <2e-16 ***
## delay_cc1 -0.34260 0.03240 -10.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## delay_cc1 0.043summary(m2_delay <- glmer(correct ~ delay_cc + (delay_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc + (delay_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 47084.2 47127.5 -23537.1 47074.2 41947
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2244 -1.1013 0.4908 0.6573 1.0211
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.42730 0.6537
## delay_cc1 0.01198 0.1095 0.44
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.13294 0.05747 19.71 <2e-16 ***
## delay_cc1 -0.30347 0.03024 -10.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## delay_cc1 -0.001Now moderated by condition
#Singular for both
# summary(m1_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc*cond_cc|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# summary(m2_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc*cond_cc|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Simplified
summary(m1_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc + cond_cc|ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 43090.3 43175.7 -21535.1 43070.3 37990
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0844 -1.0664 0.4738 0.6770 1.1528
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.35727 0.5977
## delay_cc1 0.02429 0.1559 0.37
## cond_cc1 0.50501 0.7106 -0.22 -0.17
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.10896 0.05555 19.964 <2e-16 ***
## delay_cc1 -0.35536 0.03271 -10.866 <2e-16 ***
## cond_cc1 -0.17948 0.07025 -2.555 0.0106 *
## delay_cc1:cond_cc1 -0.01904 0.05743 -0.332 0.7402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1
## delay_cc1 0.019
## cond_cc1 -0.196 -0.055
## dly_cc1:c_1 0.006 -0.041 -0.226summary(m2_delay_cond <- glmer(correct ~ delay_cc*cond_cc + (delay_cc + cond_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + (delay_cc + cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 46388.5 46474.9 -23184.2 46368.5 41942
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.9100 -1.0345 0.4587 0.6672 1.1512
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.45639 0.6756
## delay_cc1 0.01049 0.1024 0.55
## cond_cc1 0.56015 0.7484 -0.09 -0.21
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.17084 0.05937 19.723 <2e-16 ***
## delay_cc1 -0.30464 0.03033 -10.043 <2e-16 ***
## cond_cc1 -0.13973 0.07013 -1.992 0.0463 *
## delay_cc1:cond_cc1 -0.13156 0.05532 -2.378 0.0174 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1
## delay_cc1 0.023
## cond_cc1 -0.085 -0.054
## dly_cc1:c_1 0.001 -0.025 -0.215Model including appropriate covariates at the highest complexity that would fit (removed trial-within-cond random slope and delay*cond slope due to singular errors in studies 2 and 1 respectively)
# summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc + cond_cc + scale(trial_within_condition)|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
#
# # Singular
# summary(m2_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc + cond_cc + scale(trial_within_condition)|ID),
# data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
# Also singular
# summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
# (delay_cc*cond_cc + 1|ID),
# data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))
summary(m1_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
(delay_cc + cond_cc |ID),
data=s1_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +
## (delay_cc + cond_cc | ID)
## Data: s1_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 42576.5 42670.5 -21277.2 42554.5 37989
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.2409 -1.0187 0.4627 0.6726 1.4465
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.36544 0.6045
## delay_cc1 0.02941 0.1715 0.36
## cond_cc1 0.51948 0.7208 -0.21 -0.18
## Number of obs: 38000, groups: ID, 125
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.11508 0.05618 19.849 <2e-16 ***
## delay_cc1 -0.36388 0.03352 -10.857 <2e-16 ***
## cond_cc1 -0.18364 0.07119 -2.579 0.0099 **
## scale(trial_within_condition) 0.27402 0.01215 22.553 <2e-16 ***
## delay_cc1:cond_cc1 -0.01022 0.05790 -0.177 0.8599
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1 0.029
## cond_cc1 -0.194 -0.064
## scl(trl_w_) 0.018 -0.019 -0.004
## dly_cc1:c_1 0.007 -0.041 -0.225 0.006summary(m2_full_delay_cond <- glmer(correct ~ delay_cc*cond_cc + scale(trial_within_condition) +
(delay_cc + cond_cc|ID),
data=s2_train, family="binomial", control = glmerControl(optimizer = "bobyqa")))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: correct ~ delay_cc * cond_cc + scale(trial_within_condition) +
## (delay_cc + cond_cc | ID)
## Data: s2_train
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik deviance df.resid
## 45878.5 45973.6 -22928.2 45856.5 41941
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.8197 -0.9779 0.4437 0.6526 1.4299
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.46804 0.6841
## delay_cc1 0.01325 0.1151 0.47
## cond_cc1 0.57474 0.7581 -0.09 -0.18
## Number of obs: 41952, groups: ID, 138
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.17603 0.06016 19.549 <2e-16 ***
## delay_cc1 -0.30940 0.03083 -10.036 <2e-16 ***
## cond_cc1 -0.14541 0.07100 -2.048 0.0405 *
## scale(trial_within_condition) 0.26368 0.01172 22.489 <2e-16 ***
## delay_cc1:cond_cc1 -0.11653 0.05571 -2.092 0.0365 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dly_c1 cnd_c1 sc(__)
## delay_cc1 0.016
## cond_cc1 -0.084 -0.052
## scl(trl_w_) 0.016 -0.013 -0.004
## dly_cc1:c_1 0.001 -0.025 -0.214 0.010car::vif(m1_full_delay_cond)## delay_cc cond_cc
## 1.007765 1.059549
## scale(trial_within_condition) delay_cc:cond_cc
## 1.000389 1.056896car::vif(m2_full_delay_cond)## delay_cc cond_cc
## 1.004262 1.051868
## scale(trial_within_condition) delay_cc:cond_cc
## 1.000268 1.049695Heterogeneity in condition effect
Calculate some bxal differences
Study 1
s1_perf <- s1_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s1 <- as.numeric(unlist(s1_perf[s1_perf$condition=="overt","m"]-
s1_perf[s1_perf$condition=="cognitive", "m"]))
# .. and test subject-level differences in performance
s1_perf_test <- s1_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s1 <- as.numeric(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]-
s1_perf_test[s1_perf_test$condition=="cognitive", "m"]))
s1_perf_overall <- s1_train %>% group_by(ID) %>% summarize(m=mean(correct))
s1_test_perf_overall <- s1_sit %>% group_by(ID) %>% summarize(m=mean(correct))
# Modest correlations across conditions
cor.test(unlist(s1_perf[s1_perf$condition=="overt","m"]), unlist(s1_perf[s1_perf$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s1_perf[s1_perf$condition == "overt", "m"]) and unlist(s1_perf[s1_perf$condition == "cognitive", "m"])
## t = 5.461, df = 123, p-value = 2.507e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2885290 0.5729171
## sample estimates:
## cor
## 0.4417538cor.test(unlist(s1_perf_test[s1_perf_test$condition=="overt","m"]), unlist(s1_perf_test[s1_perf_test$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s1_perf_test[s1_perf_test$condition == "overt", "m"]) and unlist(s1_perf_test[s1_perf_test$condition == "cognitive", "m"])
## t = 5.6112, df = 123, p-value = 1.263e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2995928 0.5809966
## sample estimates:
## cor
## 0.4514492Examining if between-condition diffs are more common at the highest level of perf, because at lower ends more noise in bx. But not much relationship
cor.test(s1_perf_overall$m, o_min_c_perf_s1)##
## Pearson's product-moment correlation
##
## data: s1_perf_overall$m and o_min_c_perf_s1
## t = 1.0074, df = 123, p-value = 0.3157
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.08652362 0.26190545
## sample estimates:
## cor
## 0.09045834ComparePars(s1_perf_overall$m, o_min_c_perf_s1, use_identity_line = 0)
Study 2
s2_perf <- s2_train %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_perf_s2 <- as.numeric(unlist(s2_perf[s2_perf$condition=="overt","m"]-
s2_perf[s2_perf$condition=="cognitive", "m"]))
s2_perf_test <- s2_sit %>% group_by(condition, ID) %>% summarize(m=mean(correct))## `summarise()` has grouped output by 'condition'. You can override using the
## `.groups` argument.o_min_c_testperf_s2 <- as.numeric(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]-
s2_perf_test[s2_perf_test$condition=="cognitive", "m"]))
s2_perf_overall <- s2_train %>% group_by(ID) %>% summarize(m=mean(correct))
s2_test_perf_overall <- s2_sit %>% group_by(ID) %>% summarize(m=mean(correct))
cor.test(unlist(s2_perf[s2_perf$condition=="overt","m"]), unlist(s2_perf[s2_perf$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s2_perf[s2_perf$condition == "overt", "m"]) and unlist(s2_perf[s2_perf$condition == "cognitive", "m"])
## t = 6.0926, df = 136, p-value = 1.074e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3207662 0.5848974
## sample estimates:
## cor
## 0.4630508cor.test(unlist(s2_perf_test[s2_perf_test$condition=="overt","m"]), unlist(s2_perf_test[s2_perf_test$condition=="cognitive","m"]))##
## Pearson's product-moment correlation
##
## data: unlist(s2_perf_test[s2_perf_test$condition == "overt", "m"]) and unlist(s2_perf_test[s2_perf_test$condition == "cognitive", "m"])
## t = 4.8527, df = 136, p-value = 3.29e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2319669 0.5180282
## sample estimates:
## cor
## 0.3841799cor.test(s2_perf_overall$m, o_min_c_perf_s2)##
## Pearson's product-moment correlation
##
## data: s2_perf_overall$m and o_min_c_perf_s2
## t = -0.0051584, df = 136, p-value = 0.9959
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1675348 0.1666748
## sample estimates:
## cor
## -0.0004423248ComparePars(s2_perf_overall$m, o_min_c_perf_s2, use_identity_line = 0)
Invalid relationships to perf diffs
s1_o_inv <- s1_train %>% group_by(ID) %>% summarize(o_inv=mean(overt_invalid))
s1_c_inv <- s1_train %>% group_by(ID) %>% summarize(c_inv=mean(cognitive_invalid))
assert(all(s1_o_inv$ID==s1_c_inv$ID))
s1_inv_diff <- s1_c_inv$c_inv-s1_o_inv$o_inv
cor.test(o_min_c_perf_s1, s1_inv_diff)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and s1_inv_diff
## t = 0.029401, df = 123, p-value = 0.9766
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1730371 0.1781756
## sample estimates:
## cor
## 0.002650977ComparePars(o_min_c_perf_s1, s1_inv_diff, use_identity_line = 0)
s2_o_inv <- s2_train %>% group_by(ID) %>% summarize(o_inv=mean(prop_invalid_overt))
s2_c_inv <- s2_train %>% group_by(ID) %>% summarize(c_inv=mean(prop_invalid_cognitive))
assert(all(s2_o_inv$ID==s2_c_inv$ID))
s2_inv_diff <- s2_c_inv$c_inv-s2_o_inv$o_inv
cor.test(o_min_c_perf_s2, s2_inv_diff)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and s2_inv_diff
## t = 2.7091, df = 136, p-value = 0.007615
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.0614962 0.3790481
## sample estimates:
## cor
## 0.2262758ComparePars(s2_inv_diff, o_min_c_perf_s2, use_identity_line = 0, xchar="cog minus overt invalid", ychar="overt minus cog correct")
Percent of pts who actually showing better performance in cognitive
Training and test study 1
# Study 1
length(which(o_min_c_perf_s1 < 0))/length(o_min_c_perf_s1)## [1] 0.376length(which(o_min_c_perf_s1 == 0))/length(o_min_c_perf_s1) # 2.4% show no diff ## [1] 0.024length(which(o_min_c_testperf_s1 < 0))/length(o_min_c_testperf_s1)## [1] 0.328length(which(o_min_c_testperf_s1 == 0))/length(o_min_c_testperf_s1) ## [1] 0.104Training and test study 2
length(which(o_min_c_perf_s2 < 0))/length(o_min_c_perf_s2)## [1] 0.3913043length(which(o_min_c_perf_s2 == 0))/length(o_min_c_perf_s2) # < 1% show no diff## [1] 0.007246377length(which(o_min_c_testperf_s2 < 0))/length(o_min_c_testperf_s2)## [1] 0.3550725length(which(o_min_c_testperf_s2 == 0))/length(o_min_c_testperf_s2) ## [1] 0.1884058Modeling Results (Note: model numbers match those from s.R but not those used in paper…)
… because the paper numbering follows the logic of how presented there, and because we culled poor/uninformative/auxiliary models as developing the streamlined final set
Model 1 — Q learner with different decay (phi). Includes eps, CK, block explor parameter
Numbers from before bug fix on the use of prior (3/23/23) dont use: 58986 | 52624 | 57510 | 74680
m1_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior17352.csv")
m1_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior58319.csv")
m1_study1_eb <- rbind(m1_study1_eb_v1, m1_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m1_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior12123.csv")
m1_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior74336.csv")
m1_study2_eb <- rbind(m1_study2_eb_v1, m1_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m1_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")
# write.csv(m1_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPrior_merged.csv")cat("\nBeta Median \n")##
## Beta Medianmedian(m1_study1_eb$beta)## [1] 2.330957cat("SD \n")## SDsd(m1_study1_eb$beta)## [1] 5.393476cat("\nRL LR \n")##
## RL LRmedian(m1_study1_eb$q_LR)## [1] 0.5916992cat("SD \n")## SDsd(m1_study1_eb$q_LR)## [1] 0.2752539cat("\nphi cog \n")##
## phi cogmedian(m1_study1_eb$cog_phi)## [1] 0.1593268cat("SD \n")## SDsd(m1_study1_eb$cog_phi)## [1] 0.2250586cat("\nphi overt \n")##
## phi overtmedian(m1_study1_eb$overt_phi)## [1] 0.1328623cat("SD \n")## SDsd(m1_study1_eb$overt_phi)## [1] 0.169769cat("\nepsilon \n")##
## epsilonmedian(m1_study1_eb$epsilon)## [1] 0.005554817cat("SD \n")## SDsd(m1_study1_eb$epsilon)## [1] 0.02699492cat("\nChoice LR \n")##
## Choice LRmedian(m1_study1_eb$choice_LR)## [1] 0.1198229cat("SD \n")## SDsd(m1_study1_eb$choice_LR)## [1] 0.1222895cat("\nES \n")##
## ESmedian(m1_study1_eb$explor_scalar)## [1] 0.3783385cat("SD \n")## SDsd(m1_study1_eb$explor_scalar)## [1] 0.2089571cat("\nBeta Median \n")##
## Beta Medianmedian(m1_study2_eb$beta)## [1] 2.692798cat("SD \n")## SDsd(m1_study2_eb$beta)## [1] 4.910113cat("\nRL LR \n")##
## RL LRmedian(m1_study2_eb$q_LR)## [1] 0.5543408cat("SD \n")## SDsd(m1_study2_eb$q_LR)## [1] 0.2806cat("\nphi cog \n")##
## phi cogmedian(m1_study2_eb$cog_phi)## [1] 0.2033696cat("SD \n")## SDsd(m1_study2_eb$cog_phi)## [1] 0.2858507cat("\nphi overt \n")##
## phi overtmedian(m1_study2_eb$overt_phi)## [1] 0.157541cat("SD \n")## SDsd(m1_study2_eb$overt_phi)## [1] 0.2223507cat("\nepsilon \n")##
## epsilonmedian(m1_study2_eb$epsilon)## [1] 7.844961e-07cat("SD \n")## SDsd(m1_study2_eb$epsilon)## [1] 0.02668559cat("\nChoice LR \n")##
## Choice LRmedian(m1_study2_eb$choice_LR)## [1] 0.09913484cat("SD \n")## SDsd(m1_study2_eb$choice_LR)## [1] 0.128331cat("\nES \n")##
## ESmedian(m1_study2_eb$explor_scalar)## [1] 0.2382697cat("SD \n")## SDsd(m1_study2_eb$explor_scalar)## [1] 0.1951073ComparePars(m1_study1_eb$cog_phi, m1_study1_eb$overt_phi,
"Phi", "Cog", "Overt")
ComparePars(m1_study2_eb$cog_phi, m1_study2_eb$overt_phi,
"Phi", "Cog", "Overt")
Relationships at both train and test in studies 1 and 2
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m1_study1_eb$ID)
o_min_c_phi_s1_m1 <- m1_study1_eb$overt_phi - m1_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m1
## t = -6.9629, df = 123, p-value = 1.765e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6469195 -0.3927853
## sample estimates:
## cor
## -0.5317171cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m1
## t = -6.0126, df = 123, p-value = 1.926e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6018399 -0.3284889
## sample estimates:
## cor
## -0.4766035assert(s2_perf[s2_perf$condition=="overt", "ID"]==m1_study2_eb$ID)
o_min_c_phi_s2_m1 <- m1_study2_eb$overt_phi - m1_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m1
## t = -4.7721, df = 136, p-value = 4.642e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5133414 -0.2259169
## sample estimates:
## cor
## -0.3787242cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m1)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m1
## t = -7.1401, df = 136, p-value = 5.093e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6339502 -0.3889985
## sample estimates:
## cor
## -0.5221611Descriptives and chi-square tests
Higher median decay in cog in both
# hist(sqrt(m1_study1_eb$cog_phi), breaks=100)
# hist(m1_study1_eb$cog_phi, breaks=100)
# hist(m1_study1_eb$overt_phi, breaks=100)
cat("\n Study 1 cog\n")##
## Study 1 cogmedian(m1_study1_eb$cog_phi)## [1] 0.1593268cat("\n Study 1 overt\n")##
## Study 1 overtmedian(m1_study1_eb$overt_phi)## [1] 0.1328623m1_s1_phi <- rbind(
data.frame("phi"=m1_study1_eb$overt_phi, "cond"="overt"),
data.frame("phi"=m1_study1_eb$cog_phi, "cond"="cognitive")
)
m1_s2_phi <- rbind(
data.frame("phi"=m1_study2_eb$overt_phi, "cond"="overt"),
data.frame("phi"=m1_study2_eb$cog_phi, "cond"="cognitive")
)
m1_s1_phi$Experiment <- "Experiment 1"
m1_s2_phi$Experiment <- "Experiment 2"
m1_phi_ID_comb <- rbind(m1_s1_phi, m1_s2_phi)
phi_medians <- rbind(
data.frame(m1_s1_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 1"),
data.frame(m1_s2_phi %>% group_by(cond) %>% summarize(m=median(phi)), "Experiment"="Experiment 2")
)Difference higher for means but not as meaningful bc of skew, so plotting median
m1_s1_phi %>% group_by(cond) %>% summarize(m=mean(phi))## # A tibble: 2 × 2
## cond m
## <chr> <dbl>
## 1 cognitive 0.235
## 2 overt 0.176m1_s2_phi %>% group_by(cond) %>% summarize(m=mean(phi))## # A tibble: 2 × 2
## cond m
## <chr> <dbl>
## 1 cognitive 0.310
## 2 overt 0.209phi_plot <-ggplot(phi_medians, aes(x=cond, y=m, color=cond, group=cond)) +
geom_jitter(data = m1_phi_ID_comb,
aes(x=cond, y=phi, group=cond),
color="black", fill="gray57", size=2, pch=21, width=.15) +
geom_bar(stat="identity", fill="white", size=3, alpha=.2) +
facet_wrap(~ Experiment) + ga + ap + ft + lp + ylab(TeX('$\\phi') ) + xlab("") +
scale_color_manual(values=c("purple", "orange")) + ylim(-.03, 1.1) + tol
phi_plot
#ggsave("../paper/figs/fig_parts/phi_plot.png", phi_plot, height=5, width=10, dpi=700)Percent higher phi in studies 1 and 2
c_min_o_phi_s1_m1 <- m1_study1_eb$cog_phi - m1_study1_eb$overt_phi
c_min_o_phi_s2_m1 <- m1_study2_eb$cog_phi - m1_study2_eb$overt_phitable((c_min_o_phi_s1_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s1_m1 > 0)*1))## 1
## 0.576table((c_min_o_phi_s2_m1 > 0)*1)[2]/sum(table((c_min_o_phi_s2_m1 > 0)*1))## 1
## 0.6304348chisq.test(table((c_min_o_phi_s1_m1 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((c_min_o_phi_s1_m1 > 0) * 1)
## X-squared = 2.888, df = 1, p-value = 0.08924chisq.test(table((c_min_o_phi_s2_m1 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((c_min_o_phi_s2_m1 > 0) * 1)
## X-squared = 9.3913, df = 1, p-value = 0.00218m1_s1_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s1, "test_oc_diff"=o_min_c_testperf_s1, "phi_co_diff"=c_min_o_phi_s1_m1)
m1_s2_plot_df <- data.frame("train_oc_diff"=o_min_c_perf_s2, "test_oc_diff"=o_min_c_testperf_s2, "phi_co_diff"=c_min_o_phi_s2_m1)Training plots
a <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="black") +
geom_point(pch=21, fill="gray57", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("Proportion correct: \nOvert minus cognitive") +
xlab(TeX("")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
a## `geom_smooth()` using formula = 'y ~ x'
b <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=train_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="black") +
geom_point(pch=21, fill="gray57", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("") +
xlab(TeX("")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
b## `geom_smooth()` using formula = 'y ~ x'
Test plots
c <- ggplot(m1_s1_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
geom_point(pch=21, fill="gray79", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("Proportion correct: \nOvert minus cognitive") +
xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)
d <- ggplot(m1_s2_plot_df, aes(x=phi_co_diff, y=test_oc_diff)) +
geom_smooth(method="lm", se=FALSE, size=3, color="gray57") +
geom_point(pch=21, fill="gray79", size=6) +
stat_cor(method="pearson", size=6, label.y=.5) +
ylab("") +
xlab(TeX("$\\phi^{Cog}-\\phi^{Overt}")) +
ga + ap +
xlim(-.9, 1.05) + ylim(-.7, .65)all_vs <- (a + b) / (c + d)all_vs## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
#ggsave("../paper/figs/fig_parts/phi_vs_propdiff.png", all_vs, height=8, width=14, dpi=700)Correlation with decay effect
cor.test(m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi, ranef(m1_full_delay_cond)$ID$delay_cc1)##
## Pearson's product-moment correlation
##
## data: m1_study1_eb_v1$cog_phi + m1_study1_eb_v1$overt_phi and ranef(m1_full_delay_cond)$ID$delay_cc1
## t = -3.7766, df = 123, p-value = 0.0002463
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4712759 -0.1555427
## sample estimates:
## cor
## -0.3223455cor.test(m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi, ranef(m2_full_delay_cond)$ID$delay_cc1)##
## Pearson's product-moment correlation
##
## data: m1_study2_eb_v1$cog_phi + m1_study2_eb_v1$overt_phi and ranef(m2_full_delay_cond)$ID$delay_cc1
## t = -5.6966, df = 136, p-value = 7.244e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5646095 -0.2933251
## sample estimates:
## cor
## -0.438916ComparePars((m1_study1_eb_v1$cog_phi+m1_study1_eb_v1$overt_phi), ranef(m1_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)
ComparePars((m1_study2_eb_v1$cog_phi+m1_study2_eb_v1$overt_phi), ranef(m2_full_delay_cond)$ID$delay_cc1, use_identity_line = 0)
Model validation, main text parts
6/5/23 — reran sims to make sure completely updated/bug fixed version of par results used
s1_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__train_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_train_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__train_RunMQLearnerDiffDecayToPessPrior28414.csv")
s1_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_1_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior57088.csv")
s2_test_sim_m1_eb <-
rm(sp,
"SIM_EMPIRICAL_BAYES_study_2_train_SIT__sit_RunMQLearnerDiffDecayToPessPrior28414.csv")a <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s1_train, s1_train_sim_m1_eb))
b <- AltOneSimEmpPlotRep(AltPlotSimTrainingCurves(emp_df=s2_train, s2_train_sim_m1_eb))c <- PlotSimEmpTest(emp_df=s1_sit, sim_df=s1_test_sim_m1_eb)
d <- PlotSimEmpTest(emp_df=s2_sit, sim_df=s2_test_sim_m1_eb)a
b
sim_test_comb <- c + d #(temporarily flipped tol -> lp to get the legend for paper)
sim_test_comb
#ggsave("../paper/figs/fig_parts/train_sim_experiment1_plot.png", a, height=7, width=12, dpi=700)ggsave("../paper/figs/fig_parts/test_sim_experiment1_plot.png", c, height=5, width=11, dpi=700)ggsave("../paper/figs/fig_parts/test_sim_experiment2_plot.png", d, height=5, width=11, dpi=700)For supplemental because of space constraints
#ggsave("../paper/figs/fig_parts/train_sim_experiment2_plot.png", b, height=7, width=12, dpi=700)Model 3 — Same as m1 but not letting decay vary
Prior to vactoin — bug fixed files swapped in (see to_do notes)
m3_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior35433.csv")
m3_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior22212.csv")
m3_study1_eb <- rbind(m3_study1_eb_v1, m3_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m3_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior29818.csv")
m3_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior56742.csv")
m3_study2_eb <- rbind(m3_study2_eb_v1,m3_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m3_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")
# write.csv(m3_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPrior_merged.csv")Worse in AIC
ComparePars(m3_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m3", "m1")
ComparePars(m3_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m3", "m1")
sum(m3_study1_eb$AIC-m1_study1_eb$AIC)## [1] 134.026sum(m3_study2_eb$AIC-m1_study2_eb$AIC)## [1] 147.2234Sanity check that the amount of model improvement when letting decay varies correlates with the diff in phi between conditions
assert(m1_study1_eb$ID==m3_study1_eb$ID)
assert(m1_study2_eb$ID==m3_study2_eb$ID)cor.test(m3_study1_eb$AIC-m1_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), method="spearman")##
## Spearman's rank correlation rho
##
## data: m3_study1_eb$AIC - m1_study1_eb$AIC and abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi)
## S = 107038, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6711582cor.test(m3_study2_eb$AIC-m1_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), method="spearman")##
## Spearman's rank correlation rho
##
## data: m3_study2_eb$AIC - m1_study2_eb$AIC and abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi)
## S = 163526, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6266436(Function prints Pearson’s r but more appropriate Spearman’s \(\rho\) reported one cell above)
ComparePars(m1_study1_eb$AIC-m3_study1_eb$AIC, abs(m1_study1_eb$overt_phi - m1_study1_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)
ComparePars(m1_study2_eb$AIC-m3_study2_eb$AIC, abs(m1_study2_eb$overt_phi - m1_study2_eb$cog_phi), "", "AIC change", "Overt min Cog phi", use_identity_line = 0)
Model 4 — Same as m3 but to 0 inits
m4_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits47402.csv")
m4_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits84391.csv")
m4_study1_eb <- rbind(m4_study1_eb_v1, m4_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m4_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits41957.csv")
m4_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits37660.csv")
m4_study2_eb <- rbind(m4_study2_eb_v1,
m4_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m4_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")
# write.csv(m4_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayTo0Inits_merged.csv")Worse in AIC than pessimistic priors
ComparePars(m4_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m3", "m1")
ComparePars(m4_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m3", "m1")
sum(m4_study1_eb$AIC-m3_study1_eb$AIC)## [1] 120.8681sum(m4_study2_eb$AIC-m3_study2_eb$AIC)## [1] 446.3359Model 5 — Same as m1 but lets beta vary instead of decay
4/6 — bug fixed
m5_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta67330.csv")
m5_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta37484.csv")
m5_study1_eb <- rbind(m5_study1_eb_v1, m5_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m5_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta76602.csv")
m5_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffBeta11129.csv")
m5_study2_eb <- rbind(m5_study2_eb_v1, m5_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Worse in AIC than diff decay
ComparePars(m5_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m5", "m1")
ComparePars(m5_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m5", "m1")
sum(m5_study1_eb$AIC-m1_study1_eb$AIC)## [1] 672.2797sum(m5_study2_eb$AIC-m1_study2_eb$AIC)## [1] 857.6129Model 6 — Same as m1 but lets learning rate vary instead of decay
4/6 — bug fix files swapped in
m6_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR12095.csv")
m6_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR13604.csv")
m6_study1_eb <- rbind(m6_study1_eb_v1, m6_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m6_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR69123.csv")
m6_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDecayToPessPriorDiffLR43205.csv")
m6_study2_eb <- rbind(m6_study2_eb_v1, m6_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Worse in AIC than decay allowed to vary by condition
ComparePars(m6_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m6", "m1")
ComparePars(m6_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m6", "m1")
sum(m6_study1_eb$AIC-m1_study1_eb$AIC)## [1] 131.9173sum(m6_study2_eb$AIC-m1_study2_eb$AIC)## [1] 173.9311Higher overt learning rate in both
median(m6_study1_eb$cog_q_LR)## [1] 0.5048065median(m6_study1_eb$overt_q_LR)## [1] 0.556566median(m6_study2_eb$cog_q_LR)## [1] 0.4577678median(m6_study1_eb$overt_q_LR)## [1] 0.556566Learning rate correlated across conditions
ComparePars(m6_study1_eb$cog_q_LR, m6_study1_eb$overt_q_LR,
"q_LR", "Cog", "Overt")
ComparePars(m6_study2_eb$cog_q_LR, m6_study2_eb$overt_q_LR,
"q_LR", "Cog", "Overt")
Correlated with perf diffs but not as strongly as phi
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m6_study1_eb$ID)
o_min_c_q_LR_s1_m6 <- m6_study1_eb$overt_q_LR - m6_study1_eb$cog_q_LR
cor.test(o_min_c_perf_s1, o_min_c_q_LR_s1_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_q_LR_s1_m6
## t = 4.361, df = 123, p-value = 2.708e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2034049 0.5088491
## sample estimates:
## cor
## 0.3659413cor.test(o_min_c_testperf_s1, o_min_c_q_LR_s1_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_q_LR_s1_m6
## t = 2.6939, df = 123, p-value = 0.008049
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06303996 0.39525878
## sample estimates:
## cor
## 0.2360345assert(s2_perf[s2_perf$condition=="overt", "ID"]==m6_study2_eb$ID)
o_min_c_q_LR_s2_m6 <- m6_study2_eb$overt_q_LR - m6_study2_eb$cog_q_LR
cor.test(o_min_c_perf_s2, o_min_c_q_LR_s2_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_q_LR_s2_m6
## t = 3.83, df = 136, p-value = 0.0001949
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1528938 0.4553870
## sample estimates:
## cor
## 0.3120266cor.test(o_min_c_testperf_s2, o_min_c_q_LR_s2_m6)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_q_LR_s2_m6
## t = 5.4029, df = 136, p-value = 2.849e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2724024 0.5489175
## sample estimates:
## cor
## 0.420372chisq.test(table((o_min_c_q_LR_s1_m6 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((o_min_c_q_LR_s1_m6 > 0) * 1)
## X-squared = 6.728, df = 1, p-value = 0.009491chisq.test(table((o_min_c_q_LR_s2_m6 > 0)*1))##
## Chi-squared test for given probabilities
##
## data: table((o_min_c_q_LR_s2_m6 > 0) * 1)
## X-squared = 1.8551, df = 1, p-value = 0.1732Model 11 — Simple q-learner model
m11_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner24681.csv")
m11_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearner77403.csv")
m11_study1_eb <- rbind(m11_study1_eb_v1,m11_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m11_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner18790.csv")
m11_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearner26119.csv")
m11_study2_eb <- rbind(m11_study2_eb_v1, m11_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))Much worse in AIC to not include decay
ComparePars(m11_study1_eb$AIC, m3_study1_eb$AIC, "AIC", "m11", "m3")
ComparePars(m11_study2_eb$AIC, m3_study2_eb$AIC, "AIC", "m11", "m3")
sum(m11_study1_eb$AIC-m3_study1_eb$AIC)## [1] 918.2861sum(m11_study2_eb$AIC-m3_study2_eb$AIC)## [1] 1225.945Model 12 — Simple Bayes model
m12_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner13415.csv")
m12_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner42048.csv")
m12_study1_eb <- rbind(m12_study1_eb_v1,m12_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m12_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner22977.csv")
m12_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearner23536.csv")
m12_study2_eb <- rbind(m12_study2_eb_v1, m12_study2_eb_v2) %>% group_by(ID) %>% slice(which.min(nll))In turn much worse than the simple Q-learning model
ComparePars(m12_study1_eb$AIC, m11_study1_eb$AIC, "AIC", "m12", "m11")
ComparePars(m12_study2_eb$AIC, m11_study2_eb$AIC, "AIC", "m12", "m11")
sum(m12_study1_eb$AIC-m11_study1_eb$AIC)## [1] 622.5435sum(m12_study2_eb$AIC-m11_study2_eb$AIC)## [1] 1012.796Model 13 — M12 but allowing beta to vary
m13_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta60557.csv")
m13_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta85855.csv")
m13_study1_eb <- rbind(m13_study1_eb_v1, m13_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m13_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMBayesLearnerDiffBeta29859.csv")
m13_study2_eb <- rbind(m13_study2_eb_v1, m13_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Way worse than m1
ComparePars(m13_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m13", "m1")
ComparePars(m13_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m13", "m1")
sum(m13_study1_eb$AIC-m1_study1_eb$AIC)## [1] 2081.633sum(m13_study2_eb$AIC-m1_study2_eb$AIC)## [1] 2997.104Higher median beta in overt
median(m13_study1_eb$cog_beta)## [1] 1.790819median(m13_study1_eb$overt_beta)## [1] 1.876175median(m13_study2_eb$cog_beta)## [1] 1.893116median(m13_study2_eb$overt_beta)## [1] 2.099269Modest correlations between beta and perf
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m13_study1_eb$ID)
o_min_c_beta_s1_m13 <- m13_study1_eb$overt_beta - m13_study1_eb$cog_beta
cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_beta_s1_m13
## t = 3.5966, df = 123, p-value = 0.0004653
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1404791 0.4592086
## sample estimates:
## cor
## 0.3084768cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_beta_s1_m13
## t = 2.7331, df = 123, p-value = 0.007199
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06646281 0.39815509
## sample estimates:
## cor
## 0.2392776assert(s2_perf[s2_perf$condition=="overt", "ID"]==m13_study2_eb$ID)
o_min_c_beta_s2_m13 <- m13_study2_eb$overt_beta - m13_study2_eb$cog_beta
cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_beta_s2_m13
## t = 2.8955, df = 136, p-value = 0.004412
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.07696563 0.39227998
## sample estimates:
## cor
## 0.2409713cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m13)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_beta_s2_m13
## t = 1.9265, df = 136, p-value = 0.05613
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.004233115 0.321339865
## sample estimates:
## cor
## 0.162987Model 14 — Hybrid model — integrative plus bayes — bayes contribution varies by cond
Alt refers to this being the simpler version with uncertainty factored in at the action (arm) level rather than a weighted uncertainty computed — this and the more complex version were v similar in terms of fit
m14_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt39650.csv")
m14_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt71571.csv")
m14_study1_eb <- rbind(m14_study1_eb_v1, m14_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m14_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt38226.csv")
m14_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayesAlt28286.csv")
m14_study2_eb <- rbind(m14_study2_eb_v1, m14_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Original version with the more complicated weighting scheme
m14_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes25885.csv")
m14_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes65430.csv")
m14_study1_eb_old <- rbind(m14_study1_eb_v1_old, m14_study1_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))
m14_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes30853.csv")
m14_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDecayingQLearnerDiffBayes58630.csv")
m14_study2_eb_old <- rbind(m14_study2_eb_v1_old, m14_study2_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))As expected, these are very similar, with the new version a tad better
ComparePars(m14_study1_eb$AIC, m14_study1_eb_old$AIC, "AIC", "m14", "m14 old")
sum(m14_study1_eb$AIC-m14_study1_eb_old$AIC)## [1] -18.67696ComparePars(m14_study2_eb$AIC, m14_study2_eb_old$AIC, "AIC", "m14", "m14 old")
sum(m14_study2_eb$AIC-m14_study2_eb_old$AIC)## [1] -43.19563Worse AIC (and nll) compared to m1 (not a nested model bc doesn’t have diff decay)
ComparePars(m14_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m14", "m1")
ComparePars(m14_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m14", "m1")
sum(m14_study1_eb$AIC-m1_study1_eb$AIC)## [1] 448.2412sum(m14_study2_eb$AIC-m1_study2_eb$AIC)## [1] 486.6508sum(m14_study1_eb$nll-m1_study1_eb$nll)## [1] 99.12058sum(m14_study2_eb$nll-m1_study2_eb$nll)## [1] 105.3254Higher median rho in overt
median(m14_study1_eb$rho_cog)## [1] 0.1651074median(m14_study1_eb$rho_overt)## [1] 0.2611568median(m14_study2_eb$rho_cog)## [1] 0.1066345median(m14_study2_eb$rho_overt)## [1] 0.2367347Correlations w perf but not as strongly as phi
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m14_study1_eb$ID)
o_min_c_beta_s1_m14 <- m14_study1_eb$rho_overt - m14_study1_eb$rho_cog
cor.test(o_min_c_perf_s1, o_min_c_beta_s1_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_beta_s1_m14
## t = 4.2195, df = 123, p-value = 4.71e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1919767 0.4999821
## sample estimates:
## cor
## 0.3555962cor.test(o_min_c_testperf_s1, o_min_c_beta_s1_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_beta_s1_m14
## t = 3.3007, df = 123, p-value = 0.001262
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1154251 0.4388741
## sample estimates:
## cor
## 0.2852507assert(s2_perf[s2_perf$condition=="overt", "ID"]==m14_study2_eb$ID)
o_min_c_beta_s2_m14 <- m14_study2_eb$rho_overt - m14_study2_eb$rho_cog
cor.test(o_min_c_perf_s2, o_min_c_beta_s2_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_beta_s2_m14
## t = 5.1757, df = 136, p-value = 7.977e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2558970 0.5363993
## sample estimates:
## cor
## 0.4056554cor.test(o_min_c_testperf_s2, o_min_c_beta_s2_m14)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_beta_s2_m14
## t = 3.4613, df = 136, p-value = 0.0007187
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1232954 0.4311429
## sample estimates:
## cor
## 0.2845379Model 15 — Hybrid model — integrative plus bayes — decay contribution varies by cond
Alternate version with the simpler action-wise weighting (see m14 description)
m15_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt33317.csv")
m15_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt35801.csv")
m15_study1_eb <- rbind(m15_study1_eb_v1, m15_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m15_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt64665.csv")
m15_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayesAlt54458.csv")
m15_study2_eb <- rbind(m15_study2_eb_v1, m15_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))Original version with the more complicated weighting scheme
bug numbers: 75340 81994 63931 40889
m15_study1_eb_v1_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes50404.csv")
m15_study1_eb_v2_old <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes56900.csv")
m15_study1_eb_old <- rbind(m15_study1_eb_v1_old, m15_study1_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))
m15_study2_eb_v1_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes68726.csv")
m15_study2_eb_v2_old <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMHybridDiffDecayQLearnerBayes30103.csv")
m15_study2_eb_old <- rbind(m15_study2_eb_v1_old, m15_study2_eb_v2_old) %>%
group_by(ID) %>% slice(which.min(nll))As expected, again very similar — so using new/simpler one going forward
ComparePars(m15_study1_eb$AIC, m15_study1_eb_old$AIC, "AIC", "m15", "m15 old")
sum(m15_study1_eb$AIC-m15_study1_eb_old$AIC)## [1] 11.82675ComparePars(m15_study2_eb$AIC, m15_study2_eb_old$AIC, "AIC", "m15", "m15 old")
sum(m15_study2_eb$AIC-m15_study2_eb_old$AIC)## [1] -53.60774Worse AIC compared to m1
ComparePars(m15_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m15", "m1")
ComparePars(m15_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m15", "m1")
sum(m15_study1_eb$AIC-m1_study1_eb$AIC)## [1] 169.3991sum(m15_study2_eb$AIC-m1_study2_eb$AIC)## [1] 127.4026# sum(m15_study1_eb$nll-m1_study1_eb$nll)
# sum(m15_study2_eb$nll-m1_study2_eb$nll)Higher median decay in cog in both
median(m15_study1_eb$cog_phi)## [1] 0.2204569median(m15_study1_eb$overt_phi)## [1] 0.1742236median(m15_study2_eb$cog_phi)## [1] 0.3308409median(m15_study1_eb$overt_phi)## [1] 0.1742236ComparePars(m15_study1_eb$cog_phi, m15_study1_eb$overt_phi,
"Phi", "Cog", "Overt")
ComparePars(m15_study2_eb$cog_phi, m15_study2_eb$overt_phi,
"Phi", "Cog", "Overt")
Similar perf correlations as m1
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m15_study1_eb$ID)
o_min_c_phi_s1_m15 <- m15_study1_eb$overt_phi - m15_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m15
## t = -7.1705, df = 123, p-value = 6.097e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6560039 -0.4060503
## sample estimates:
## cor
## -0.542943cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m15
## t = -6.1734, df = 123, p-value = 8.902e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6098839 -0.3397791
## sample estimates:
## cor
## -0.4863663assert(s2_perf[s2_perf$condition=="overt", "ID"]==m15_study2_eb$ID)
o_min_c_phi_s2_m15 <- m15_study2_eb$overt_phi - m15_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m15
## t = -3.8478, df = 136, p-value = 0.0001825
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4565342 -0.1543079
## sample estimates:
## cor
## -0.3133333cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m15)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m15
## t = -6.5236, df = 136, p-value = 1.247e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6058791 -0.3496199
## sample estimates:
## cor
## -0.4882024Model 27 — Fixing both ES and epsilon at the group level
m27_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed46645.csv")
m27_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed33416.csv")
m27_study1_eb <- rbind(m27_study1_eb_v1, m27_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m27_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed20328.csv")
m27_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed44601.csv")
m27_study2_eb <- rbind(m27_study2_eb_v1, m27_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m27_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")
# write.csv(m27_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorESAndEpsFixed_merged.csv")Does worsen NLL substantially
ComparePars(m27_study1_eb$nll, m1_study1_eb$nll, "nll", "m27", "m1")
ComparePars(m27_study2_eb$nll, m1_study2_eb$nll, "nll", "m27", "m1")
sum(m27_study1_eb$nll-m1_study1_eb$nll)## [1] 259.0197sum(m27_study2_eb$nll-m1_study2_eb$nll)## [1] 367.588Still higher median decay in cog in both
median(m27_study1_eb$cog_phi)## [1] 0.1637271median(m27_study1_eb$overt_phi)## [1] 0.1337206median(m27_study2_eb$cog_phi)## [1] 0.2304866median(m27_study1_eb$overt_phi)## [1] 0.1337206Strong relationships at both train and test in studies 1 and 2
assert(s1_perf[s1_perf$condition=="overt", "ID"]==m27_study1_eb$ID)
o_min_c_phi_s1_m27 <- m27_study1_eb$overt_phi - m27_study1_eb$cog_phi
cor.test(o_min_c_perf_s1, o_min_c_phi_s1_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s1 and o_min_c_phi_s1_m27
## t = -6.9726, df = 123, p-value = 1.68e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6473518 -0.3934142
## sample estimates:
## cor
## -0.5322504cor.test(o_min_c_testperf_s1, o_min_c_phi_s1_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s1 and o_min_c_phi_s1_m27
## t = -6.4826, df = 123, p-value = 1.966e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6248649 -0.3610156
## sample estimates:
## cor
## -0.504631assert(s2_perf[s2_perf$condition=="overt", "ID"]==m27_study2_eb$ID)
o_min_c_phi_s2_m27 <- m27_study2_eb$overt_phi - m27_study2_eb$cog_phi
cor.test(o_min_c_perf_s2, o_min_c_phi_s2_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_perf_s2 and o_min_c_phi_s2_m27
## t = -4.8922, df = 136, p-value = 2.775e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5203133 -0.2349247
## sample estimates:
## cor
## -0.3868433cor.test(o_min_c_testperf_s2, o_min_c_phi_s2_m27)##
## Pearson's product-moment correlation
##
## data: o_min_c_testperf_s2 and o_min_c_phi_s2_m27
## t = -7.114, df = 136, p-value = 5.847e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6328069 -0.3873771
## sample estimates:
## cor
## -0.52077159 and 62% of pts show a higher phi in overt
table((o_min_c_phi_s1_m27 > 0)*1)##
## 0 1
## 73 52table((o_min_c_phi_s2_m27 > 0)*1)##
## 0 1
## 86 52table((o_min_c_phi_s1_m27 > 0))[1]/sum(table((o_min_c_phi_s1_m27 > 0)))## FALSE
## 0.584table((o_min_c_phi_s2_m27 > 0))[1]/sum(table((o_min_c_phi_s2_m27 > 0)))## FALSE
## 0.6231884chisq.test(table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0)))##
## Chi-squared test for given probabilities
##
## data: table(if_else((o_min_c_phi_s1_m27 > 0) == TRUE, 1, 0))
## X-squared = 3.528, df = 1, p-value = 0.06034chisq.test(table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0)))##
## Chi-squared test for given probabilities
##
## data: table(if_else((o_min_c_phi_s2_m27 > 0) == TRUE, 1, 0))
## X-squared = 8.3768, df = 1, p-value = 0.0038Model 29 — Model 1 but with no choice kernel
m29_study1_eb_v1 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK21030.csv")
m29_study1_eb_v2 <- rm(bp, "BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK89190.csv")
# Replace when back
#m29_study1_eb <- m29_study1_eb_v1
m29_study1_eb <- rbind(m29_study1_eb_v1, m29_study1_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
m29_study2_eb_v1 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK64103.csv")
m29_study2_eb_v2 <- rm(bp, "BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK74167.csv")
# Delete when final one is back
m29_study2_eb <- m29_study2_eb_v1
m29_study2_eb <- rbind(m29_study2_eb_v1, m29_study2_eb_v2) %>%
group_by(ID) %>% slice(which.min(nll))
# write.csv(m29_study1_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_1_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")
# write.csv(m29_study2_eb,
# "../model_res/opts_mle_paper_final/best/BEST_study_2_train_SIT_EMPIRICAL_BAYES_RunMQLearnerDiffDecayToPessPriorNoCK_merged.csv")Dramatically worse fit
ComparePars(m29_study1_eb$AIC, m1_study1_eb$AIC, "AIC", "m29", "m1")
ComparePars(m29_study2_eb$AIC, m1_study2_eb$AIC, "AIC", "m29", "m1")
sum(m29_study1_eb$AIC-m1_study1_eb$AIC)## [1] 1815.325sum(m29_study2_eb$AIC-m1_study2_eb$AIC)## [1] 2777.038Model Comparison
AIC_diff_df <-
rbind(
data.frame("diff"=m12_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes"),
data.frame("diff"=m12_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes"),
data.frame("diff"=m4_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decayto0"),
data.frame("diff"=m4_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decayto0"),
data.frame("diff"=m3_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay"),
data.frame("diff"=m3_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay"),
data.frame("diff"=m5_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_beta"),
data.frame("diff"=m5_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_beta"),
data.frame("diff"=m13_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Bayes_diff_beta"),
data.frame("diff"=m13_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Bayes_diff_beta"),
data.frame("diff"=m6_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_lr"),
data.frame("diff"=m6_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_lr"),
data.frame("diff"=m14_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_bayes"),
data.frame("diff"=m14_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_bayes"),
data.frame("diff"=m15_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_diff_decay_bayes"),
data.frame("diff"=m15_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_diff_decay_bayes"),
data.frame("diff"=m1_study1_eb$AIC-m11_study1_eb$AIC, "study"=1, "model"="Q_decay_diff_phi"),
data.frame("diff"=m1_study2_eb$AIC-m11_study2_eb$AIC, "study"=2, "model"="Q_decay_diff_phi")
)
AIC_diff_df$study <- as.factor(AIC_diff_df$study)
AIC_diff_summ <- AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff))## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.sum(m1_study1_eb$AIC-m3_study1_eb$AIC)## [1] -134.026sum(m1_study2_eb$AIC-m3_study2_eb$AIC)## [1] -147.2234Confirm min in both studies is diff phi
AIC_diff_summ[AIC_diff_summ$study == 1, ] %>% arrange(m)## # A tibble: 9 × 3
## # Groups: model [9]
## model study m
## <chr> <fct> <dbl>
## 1 Q_decay_diff_phi 1 -8.42
## 2 Q_decay_diff_lr 1 -7.36
## 3 Q_decay 1 -7.35
## 4 Q_diff_decay_bayes 1 -7.06
## 5 Q_decayto0 1 -6.38
## 6 Q_decay_diff_bayes 1 -4.83
## 7 Q_decay_diff_beta 1 -3.04
## 8 Bayes 1 4.98
## 9 Bayes_diff_beta 1 8.23AIC_diff_summ[AIC_diff_summ$study == 2, ] %>% arrange(m)## # A tibble: 9 × 3
## # Groups: model [9]
## model study m
## <chr> <fct> <dbl>
## 1 Q_decay_diff_phi 2 -9.95
## 2 Q_diff_decay_bayes 2 -9.03
## 3 Q_decay 2 -8.88
## 4 Q_decay_diff_lr 2 -8.69
## 5 Q_decay_diff_bayes 2 -6.42
## 6 Q_decayto0 2 -5.65
## 7 Q_decay_diff_beta 2 -3.74
## 8 Bayes 2 7.34
## 9 Bayes_diff_beta 2 11.8levels <-
c(
"Bayes_diff_beta",
"Bayes",
"Q_decayto0",
"Q_decay",
"Q_decay_diff_beta",
"Q_decay_diff_lr",
"Q_decay_diff_bayes",
"Q_diff_decay_bayes",
"Q_decay_diff_phi"
)
assert(length(levels)==9)
AIC_diff_summ$model <- factor(AIC_diff_summ$model, levels = levels)
AIC_diff_df$model <- factor(AIC_diff_df$model, levels = levels)mp1 <- ggplot(AIC_diff_summ %>% filter(study == 1), aes(x=model, y=m, group=model, fill=model)) +
geom_hline(yintercept = 0, size=2) +
geom_jitter(data=AIC_diff_df %>% filter(study == 1), aes(x=model, y=diff, group=model),
fill="gray57", pch=21, width=.3, size=4, alpha=.05) +
geom_bar(stat="identity", color="black", alpha=.8, size=2) +
ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
tol + ga + ap + tp +
ggtitle("Study 1") +
#facet_wrap(~ study) +
xlab("") +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
theme(axis.title.y=element_text(size=22)) + ylim(-65, 65) + coord_flip()mp2 <- ggplot(AIC_diff_summ %>% filter(study == 2), aes(x=model, y=m, group=model, fill=model)) +
geom_hline(yintercept = 0, size=2) +
geom_jitter(data=AIC_diff_df %>% filter(study == 2), aes(x=model, y=diff, group=model),
fill="gray57", pch=21, width=.3, size=4, alpha=.05) +
geom_bar(stat="identity", color="black", alpha=.8, size=2) +
ylab(TeX('$\\Delta\ AIC$ from simple Q-learner')) +
tol + ga + ap + tp +
ggtitle("Study 2") +
xlab("") +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
theme(axis.title.y=element_text(size=22)) + ylim(-65, 65) + coord_flip()combined_aic <- mp1 + mp2combined_aic
#ggsave("../paper/figs/fig_parts/aic_plot.png", combined_aic, height=5, width=10, dpi=700)AIC_diff_df %>% group_by(model, study) %>% summarize(m=mean(diff)) %>% arrange(m)## `summarise()` has grouped output by 'model'. You can override using the
## `.groups` argument.## # A tibble: 18 × 3
## # Groups: model [9]
## model study m
## <fct> <fct> <dbl>
## 1 Q_decay_diff_phi 2 -9.95
## 2 Q_diff_decay_bayes 2 -9.03
## 3 Q_decay 2 -8.88
## 4 Q_decay_diff_lr 2 -8.69
## 5 Q_decay_diff_phi 1 -8.42
## 6 Q_decay_diff_lr 1 -7.36
## 7 Q_decay 1 -7.35
## 8 Q_diff_decay_bayes 1 -7.06
## 9 Q_decay_diff_bayes 2 -6.42
## 10 Q_decayto0 1 -6.38
## 11 Q_decayto0 2 -5.65
## 12 Q_decay_diff_bayes 1 -4.83
## 13 Q_decay_diff_beta 2 -3.74
## 14 Q_decay_diff_beta 1 -3.04
## 15 Bayes 1 4.98
## 16 Bayes 2 7.34
## 17 Bayes_diff_beta 1 8.23
## 18 Bayes_diff_beta 2 11.8Get AICs for supplemental table
cat(round(sum(m1_study1_eb$AIC), 0), "\n")## 42125cat(round(sum(m1_study2_eb$AIC), 0))## 42798cat(round(sum(m3_study1_eb$AIC), 0), "\n")## 42259cat(round(sum(m3_study2_eb$AIC), 0))## 42945cat(round(sum(m4_study1_eb$AIC), 0), "\n")## 42380cat(round(sum(m4_study2_eb$AIC), 0))## 43391cat(round(sum(m5_study1_eb$AIC), 0), "\n")## 42797cat(round(sum(m5_study2_eb$AIC), 0))## 43655cat(round(sum(m6_study1_eb$AIC), 0), "\n")## 42257cat(round(sum(m6_study2_eb$AIC), 0))## 42972cat(round(sum(m13_study1_eb$AIC), 0), "\n")## 44206cat(round(sum(m13_study2_eb$AIC), 0))## 45795cat(round(sum(m14_study1_eb$AIC), 0), "\n")## 42573cat(round(sum(m14_study2_eb$AIC), 0))## 43284cat(round(sum(m15_study1_eb$AIC), 0), "\n")## 42294cat(round(sum(m15_study2_eb$AIC), 0))## 42925cat(round(sum(m11_study1_eb$AIC), 0), "\n")## 43177cat(round(sum(m11_study2_eb$AIC), 0))## 44171cat(round(sum(m12_study1_eb$AIC), 0), "\n")## 43800cat(round(sum(m12_study2_eb$AIC), 0))## 45184